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Abstract 

Population genetic structure, historical biogeography and historical demography of the alpine toad Scutiger ningshanensis 
were studied using the combined data mtDNA cytochrome b (cyt b) and the mtDNA cytochrome c oxidase subunit I (COI) as 
the molecular markers. This species has high genetic variation. There was a significant genetic differentiation among most 
populations. Three lineages were detected. The phylogenetic relationship analyses and the SAMOVA (spatial analysis of 
molecular variance) results showed significant phylogeographic structure. 82.15% genetic variation occurred among 
populations whereas differentiation within populations only contributed 1 7.85% to the total. Mantel test results showed a 
significant correlation between the pairwise calculated genetic distance and pairwise calculated geographical distance of 
the populations (regression coefficient =0.001286, correlation coefficient =0.77051, p (r rand >r obs ) =0.0185<0.05), 
indicating the existence of isolation-by-distance pattern of genetic divergence for cyt b + COI sequence, which suggests 
that the distribution of genetic variation is due to geographical separation rather than natural selection. The population 
expansion or contraction and genetic differentiation between populations or lineages could be explained by topography 
and the repetitive uplifts of the Tsinling Mountains and the climatic cycles during the late Pliocene and Pleistocene. S. 
ningshanensis experienced a rapid population expansion about 40,000 years before present. The current decline in 
population size was probably caused by anthropogenic disturbance. Current populations of S. ningshanensis are from 
different refugia though the location of these refugia could not be determined in our study. Topography, climatic changes 
and repetitive population expansion/contraction together led to the high level of genetic variation in 5. ningshanensis. A 
total of three management units (MUs) was determined, which must be considered when conservation policy is made in the 
future. 
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Introduction 

Population genetic structure refers to the geographical pattern 
of genetic diversity within or among populations. It could be 
influenced by gene flow, genetic drift, selection, mutation and 
recombination. Gene flow is caused by the movement of 
individuals from one population to another [1]. Estimation of 
the gene flow level allows conservation biologists to understand the 
relationships between populations and assess levels of genetic 
variation in order to evaluate the relative levels of conservation 
concern hierarchically across populations in a species. Genetic 
drift is the change in the frequency of a gene variant in a 
population due to random sampling [2]. http://en.wikipedia.org/ 
wiki/Genetic_drift - cite_note-Masel_201 1 - 1 ^citejiote-Ma- 
sel_201 1-1 Genetic drift may lead to disappearance of gene 
variants and thereby reduce genetic diversity. 

Phylogeography connects historical processes in evolution with 
spatial distributions [4] . Analysis of mitochondrial data promoted 
the empirical development of phylogeography [3]. The statistical 
phylogeography is one of the widely used approaches in 
phylogeography, which takes into account the stochasticity of 



genetic processes into demographic inference based on coalescent 
models for parameter estimation [4,5]. 

The Tsinling Mountains are located in the central part of 
China, stretching from west to east (Fig. 1). These Mountains are 
boundary between Oriental realm and Palaearctic realm accord- 
ing to the zoogeographical regions of China [6], and also the 
watershed for Yangtze River and Yellow River catchment areas, 
for climate, flora and fauna in China [7-9]. The Oriental and 
Palaearctic species congregate here forming a specific biotic 
province and containing rich animal and plant resources [10,1 1]. 
Like other regions in the northern hemisphere, the Tsinling 
Mountains experienced several glacial-interglacial cycles during 
Pleistocene [12-18]. The climate associated with Pleistocene 
glacial cycles in East Asia was likely mild and characterized by a 
mosaic of mountains [19,20]. The past climatic events, such as the 
Quaternary glaciation, are believed to have played an important 
role in forming the geographical pattern of the montane species 
and could leave the vestiges in geographical distribution of genetic 
diversity of population [21-26]. The founder effect during the 
postglacial population recovery causes a reduction in population 
genetic diversity [27,28], and the subsequent rapid population 
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Figure 1. Locations of sampled populations and geographical distribution of S. ningshanensis clades on the Tsinling Mountains. Nsc 

is also the type locality of S. ningshanensis. 
doi:10.1371/journal.pone.0100729.g001 



expansion [29] may erase the previous geographical differences of 
the genetic diversity. 

The alpine toad Scutiger ningshanensis was described from the 
western part of the Tsinling Mountains [30] (Fig. 1). Four years 
later, the second specimen of this species was collected from the 
same locality [31]. Since then, other specimen of this species was 
not collected until 2009 when some specimens were collected from 
several localities in the eastern part of the Tsinling Mountains 
[32]. Other than the reports on collection of additional specimens, 
only the biological characteristics of tadpoles of this species was 
studied [33]. The habitat of this species was roughly divided into 
two parts: the western part and the eastern part. Is this 
geographical pattern caused by habitat fragmentation or by 
populations from different glacial refugia? Does the isolation by 
distance between the local populations result in occurrence of any 
speciation events? The aims of the present study were to explore 



the population genetic structure, historical biogeography and the 
historical demography of S. ningshanensis. 

Materials and Methods 

Ethics statement 

This study was approved by the Institutional Animal Care and 
Use Committee (IACUC) of Shaanxi Normal University. Only the 
clipped toes or tail tips of tadpoles were used for extraction of total 
DNA. No specific permissions were required for these locations/ 
activities. This species was ranked "Endangered B2ab(iii)" in "The 
IUCN Red List of Threatened Species™ 2013" based on an 
outdated information that this species was only found at the type 
locality (nsc) [31] (http://www.iucnredlist.org). Actually, in 
addition to the type locality, this species was reported later in 
2009 from a variety of localities including the Baiyunshan 
Mountains (bys), Shirenshan Mountains (srs) and Laojunshan 



Table 1. Sampling information and haplotypes based on cyt b and COI for 6 sampled populations of Scutiger ningshanensis. 



Population 


Location 


n 


GPS coordinates 


Elevation (m) 


Haplotypes 


hby 


Huangbaiyuan, Taibai Co., 
Shaanxi Prov. 


18 


33.8749N 107.5168E 


1652 


hbyl (1), hby10 (1), hby12 (4), hby13 (1), hby14 
(1), hby15 (1), hbyl 6 (1),hby19 (1), hby2 (1), 
hby20 (1), hby3 (1), hby4 (1), hby7 (1), hby8 (1), 
hby9 (1) 


Ify 


Liangfengya, Foping Co., 
Shaanxi Prov. 


17 


33.6668N 107.8529E 


2047 


hby13 (1), Ifyl(l), IfylO (1), Ify12(2), Ifyl 3 (1), 
Ify14 (1), Ifyl 5 (1), Ifyl 6 (1), Ifyl 7 (1), Ify2(1), Ify3 
(1), Ify4 (1), Ify5 (1), Ify6(1), Ify7 (1), Ify8(1) 


nsc 


Pingheliang, Ningshan Co., 
Shaanxi Prov. 


15 


33.4744N 108.5253E 


2000 


nscl (1), nsdO (1), nscl 1 (1), nsc12 (1), nsc13 
(1), nsc14 (1), nsc15(1), nsc2 (1), nsc3 (1), 
nsc4(1), nsc5 (1), nsc6 (1), nsc7(1), nsc8 (1), 
nsc9 (1) 


Ijs 


Laojunshan, Luanchuan Co., 
Henan Prov. 


18 


33.7272N 111.6309 


1590 


bys4 (6), Ijsl (4), IjslO (2), Ijsll (1), Ijs18 (1), Ijs4 
(1), Ijs5 (2), Ijs6 (1) 


bys 


Baiyunshan, Songxian Co., 
Henan Prov. 


16 


33.6535N 111.8283E 


1675 


bys1(8), byslO (1), bysll (2), bys14 (1), bys15 
(1), bys4 (1), bys7 (2), 


srs 


Shirenshan, Lushan Co., 
Henan Prov. 


15 


33.7286N 112.2542E 


1642 


Ijsll (5), srsl (3), srs15 (1), srs16 (1), srs17 (1), 
srs18 (1), srs5 (1), srs6 (1), srs9 (1) 



n, sample size. 

doi:1 0.1 371 /journal.pone.OI 00729.t001 
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hby-lfy nsc ljs-bys-srs 

52.7% 35.1% 12.2% 




Figure 2. The model used to test the refugial hypotheses for 5. 

ningshanensis using coalescent simulations. A single-refugium 
hypothesis concerning the refugia during the Dali glaciation (the last 
maximum glaciation in China which occurred about 50 ka before 
present) was tested. The detail interpretation for this model is given in 
the text. Branch lengths are time in generations based on a 6-year 
generation time in 5. ningshanensis. Branch widths (effective population 
size, N e ) are scaled for each group based on the proportion of the total 
N e that each group comprised. 
doi:1 0.1 371 /journal.pone.01 00729.g002 

Mountains (ljs) [32] (Fig. 1), which indicated that this species 
doesn't meet the criteria of critically endangered or endangered 
defined in The IUCN Red List of Threatened Species™ 2013, 
therefore this species should be removed from the red list. 
However, the new data was not taken into account when this 
species was ranked "Endangered B2ab(iii)" in "The IUCN Red 
List of Threatened Species™ 2013". The specific location (GPS 
coordinates) of our study was given in Table 1 . 

Sampling and laboratory protocols 

Our sampling covers the entire known distribution of this 
species. Furthermore, to make an extensive sampling, we explored 
the whole Tsinling Mountains, and fortunately collected this 
species at two locations where the distribution of this species has 
not been recorded. A total of 99 samples were collected from 6 
localities during 2011 and 2013 (Table 1, Fig. 1). Eight samples of 
the alpine toad Scutiger boukngeri were collected from Jone County 
(34.539922N 103.491647E), southern Gansu Province, China. S. 
boukngeri will be used as outgroup in phylogenetic relationship 
analysis. 

The clipped toes or tail tips of tadpoles were preserved in 100% 
ethanol and stored at — 20°C. A continuous fragment (1009 bp) of 
the mitochondrial cytochrome b (cyt b) was amplified using PCR 
(MyCycler Thermal Cycler), with primers FrogGlu-f 5'- 
TGATCTGAAAAACC ACCGTTG-3 ' and FrogThr-r 5'- 

Table 2. Genetic diversity of each population of S. ningshant 



CTCCATTCTTCGRCTTACAAG-3' [34]. A continuous frag- 
ment (63 1 bp) of the mitochondrial cytochrome c oxidase subunit 
I (COI) was amplified using PCR (MyCycler Thermal Cycler), 
with primers forward LepF5'-ATT CAA CCA ATC ATA AAG 
ATA TTG G-3' and reverse LepR5'-TAA ACT TCT GGA 
TGT CCA AAA AAT CA-3' [35]. The PCR products were 
purified using a purification kit (DC351 1-02/3514-02 250 Preps, 
Biomiga, USA). Sequencing reactions were carried out with the 
PCR primers using ABI Prism BigdyeTM Terminator Cycle 
Sequencing Ready Reaction Kit on ABI 3730XL sequencer. All 
sequences have been deposited in the GenBank databases under 
accession numbers KF757340-KF757391 (S. ningshanensis cyt b), 
KF757392-KF757439 (S. ningshanensis COI); KJ082065- 
KJ082072 (S. boukngeri cyt b), KJ082073-KJ082080 (S. boukngeri 
COI). Indicators of nuclear mitochondrial pseudogenes (numts), 
such as indels, stop codons and double peaks in sequence 
chromatograms [36] were not found. 

Determination of generation time 

A skeletochronological study of longevity of S. ningshanensis was 
conducted. The clipped phalanges were stored in 95% ethanol 
solution. The phalanges were washed in running tap water, then 
decalcified in 3% nitric acid for 12 to 24 hours. The mid- 
diaphyseal region of the phalanges was cross-sectioned at 12— 
20 u.m using a microtome and stained with hematoxylin for two 
min. Sections were examined under a light microscope and the 
number of lines of arrested growth (LAGs) was counted [37,38]. 
The number of LAGs represents the age of the toad. 

Nucleotide polymorphism 

The sequences were aligned with Clustal XI. 83 [39]. The 
aligned sequences were edited using the program BioEdit 7.0.9.0 
[40]. All analyses were performed based on the combined 
mitochondrial DNA data. Haplotype inference was conducted 
through Collapse 1.2 (http://darwin.uvigo.es). The number of 
variable and parsimony informative sites was determined using the 
program DnaSP 5.10.01 [41], and haplotype diversity (Hdj and 
nucleotide diversity (Pi) were determined through Arlequin 3.5.1.2 
[42]. 

Phylogenetic analyses 

The substitution model selection was implemented in jModelT- 
est 2.1.4 [43], the TIM2+I+G model was selected for all datasets 
by likelihood ratio tests either under the Akaike Information 
Criterion (AIC) or under the Bayesian Information Criterion 
(BIC). Bayesian inference (BI) was used to generate a phylogenetic 
hypothesis of the DNA haplotypes. BI was performed in MrBayes 





Population 


Haplotype diversity ±S.D. 


Mean number of pairwise differences ±S.D. 


Nucleotide diversity ±S.D. 


hby 


0.9591 ±0.0359 


20.1 871 35±9.337776 


0.012309±0.006363 


Ify 


0.991 7±0.0254 


14.416667±6.820745 


0.008791 ±0.004658 


nsc 


1.0000±0.0243 


8.8761 90±4.336499 


0.00541 2±0.002965 


ljs 


0.8498±0.0426 


1.928854±1.137129 


0.001 176±0.000773 


bys 


0.7500±0.1071 


2.091667±1.231968 


0.00 1275 ±0.000841 


srs 


0.9333±0.0773 


3.044444+1.729647 


0.001856±0.001193 


Total population 


0.9825 ±0.0055 


41.822511 ±18.302675 


0.025502±0.012361 



S.D., standard deviation. 

doi:1 0.1 371 /journal.pone.01 00729.t002 
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0.05 



Figure 3. Bayesian tree for the 67 sampled haplotypes of 5. 
ningshanensis based on the combined mtDNA cyt b and COI 
sequences. The Bayesian posterior probabilities from Bayesian 
analyses are presented above or under the main branches. The scale 
bar represents substitutions per site. 
doi:1 0.1 371 /journal.pone.01 00729.g003 

3. 2 [44] with 1,200,000 generations, sampling trees every 100 
generations. Two independent runs each with four simultaneous 
Monte Carlo Markov chains (MCMC) were conducted. The first 
25% of generations were discarded as 'burn-in'. The convergence 
of chains was confirmed until average standard deviation of split 
frequency is below 0.01 (0.009889) and the potential scale 



reduction factor (PSRF) is close to 1.0 for all parameters. In 
phylogenetic analysis S. boulengeri was used as outgroup. 

Population structure analyses 

The population comparisons using pairwise difference as 
distance method, and the partition of genetic diversity within 
and among populations were analyzed by analysis of molecular 
variance (AMOVA) [45] using Arlequin3.5.1.2 [42] with 10,000 
permutations. Mantel tests [46] were also conducted in Arle- 
quin3.5.1.2 to assess the significance of isolation by distance (IBD) 
between populations with 10,000 random permutations on 
matrices of pairwise population _F ST values and the geographical 
distances. Pair-wise F ST values between populations were estimat- 
ed using Arlequin3.5.1.2, while geographical distances between 
populations were calculated online at http://www.gpsvisualizer. 
com/ calculators^distance. 

The spatial genetic structure of haplotypes was analyzed using 
the program SAMOVA1.0 [47] (http://web.unife.it/progetti/ 
genetica/Isabelle/samova.html) with 1,000 permutations. The 
number of initial conditions was set to 1 00 as recommended by 
Dupanloup et al. [48] . The number K of groups of populations 
ranged from 2 to 4. The iTwifh the highest Fct represents the best 
number of groups and the best population configuration. 

Historical biogeography 

The effective population size (jV e ) of each clade (geographical 
group) for coalescent simulations was converted from Theta using 
the equation 8 = jV" e u, with u. = 0.65 x 10~ (per lineage per million 
years) x6 (generation time of S. ningshanensis). The 0-values were 
estimated using maximum likelihood method in the program 
Lamarc2.7.5 [48]. Total jV c was the sum of the jV e for all groups 
and the proportion of total jV e that each group comprised were 
used to scale the branch width of hypothesized trees (models of 
population divergence) (Fig. 2) [49-51]. Branch widths can be 
controlled by the Adjust Lineage Widths tool (the horizontal ruler) 
in the Tree Window in the program Mesquite2.75 [52]. 

A single-refugium hypothesis of population divergence during 
the Dali glaciation (the last maximum glaciation in China which 
occurred about 50 ka before present) was tested using the 
maximum likelihood estimation implemented in Mesquite2.75 
[52] to infer the distributional area of the most recent common 
ancestor (MRCA) of clades. The single-refugium model hypoth- 
esizes that all current geographical populations derived from a 
single refugium at the end of the Dali glaciation and that all 
population divergences were concurrent and resulted from the 
fragmentation of a widely distributed common ancestor's range 
(Fig. 2). Branch lengths are time in generations based on a six-year 
generation time in S. ningshanensis. The clades were treated as taxon 
states and each haplotype as character states. The haplotype trees 
by coalescence within a diverging population tree (model of 
population divergence) were simulated, and fit of a haplotype tree 
to a population tree was calculated to search for population trees 
that optimize fit of gene trees. 

Historical demography 

Neutrality tests were calculated in Arlequin3.5.1.2 [41], Fu's F & 
test [53] and Tajima's D [54] were used to detect evidence of 
recent demographic expansion within each inferred clade, under 
which negative value are expected [55] . Inference of population 
expansion events was performed using a mismatch analysis [56,57] 
using Arlequin3.5. 1.2 with the number of bootstrap replicates set 
to 10,000 to explore the demographic history of the studied 
populations. A recent growth is expected to generate a unimodal 
distribution of pairwise differences between sequences [56]. The 
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Table 3. Results of analysis of molecular 


variance 


(AMOVA) of S. ningsh 


anensis. 






Source of variation 


d.f. 


Sum of squares 


Variance components 


Percentage of variation 


Among populations 


5 


1646.756 


1 9.92646 Va 


82.15 


Within populations 


93 


402.547 


4.32847 Vb 


17.85 


Total 


98 


2049.303 


24.25492 




Fixation Index: F st = 0.82154 (p-value =0.0±0.0) 



d.f., degrees of freedom. 

doi:1 0.1 371 /joumal.pone.01 00729.t003 



validity of the expansion model was tested using the sum of 
squared deviations (SSD) and Harpending's raggedness index (R) 
between observed and expected mismatches. 

The site-frequency spectrum (for segregating sites) was calcu- 
lated in DnaSP5. 10.01 [41] to detect the excess of singleton 
mutations. The null hypothesis of the neutral model (constant 
population size) was rejected when the allelic frequency spectrum 
for the entire population revealed an excess of singleton mutations. 
The excess of singleton mutations could be caused by the 
expansion [58]. 

Furthermore, the Bayesian Skyline Plot (BSP) approach was 
used to estimate the demographic history in BEAST vl.8.0 [59]. 
The log file was analyzed in Tracer. Strict molecular clock model 
was selected based on the results from Tracer. A mean mutation 
rate of 0.65% change per lineage per million years was used [60- 
67]. Two independent BEAST runs from the same XML file were 
carried out and then the log output files were combined using 
LogCombiner. The combined log file was analyzed in Tracer to 
see if the Trace for each parameter has converged well on a 
stationary distribution. 

Detection of cryptic species 

To detect the existence of potential cryptic species, uncorrected 
p-distances [68,69] for all lineage pairs were calculated in 
PAUP4.0 [70] from all sequences. The average p-distances of all 
possible pairs of sequences (every sequence pair contains sequences 
from different lineages) were calculated. 

Results 

Generation time 

The average longevity of tadpoles was estimated by Lu et al. 
[33]. The whole longevity equals the longevity of adult plus the 
longevity of tadpole. The average duration of the tadpole stage of 
S. ningshanensis is three years [33] and the average duration of the 
adult stage after metamorphosis is six years. The average life span 



of S. ningshanensis is nine years. The average generation time (GT) 
of S. ningshanensis is six years. 

Genetic variation 

A total of 107 samples were sequenced, including 99 samples of 
S. ningshanensis and eight of S. boulengeri. A total of 1009 base pairs of 
cyt b gene and 631 base pairs of COI gene was obtained, 67 
haplotypes were identified for the combined gene sequences (cyt 
b+COI) of S. ningshanensis. Of the combined 1640 nucleotide sites, 
181 were variable (polymorphic sites or segregating sites), 121 were 
parsimony informative, and 60 were singleton variable. The 
haplotype diversity of total and most sampled population was very 
high, however, the nucleotide diversity of every sampled popula- 
tion was low (Table 2). Among the total haplotypes, 82.09% are 
private haplotypes (Table 1). 

Phylogenetic relationships among haplotypes 

BI analysis revealed three distinct clades (lineages) (hby-lfy, nsc, 
ljs-bys-srs) in S. ningshanensis. Clades hby-lfy and nsc are well 
supported with posterior probabilities of 1 and 0.99. On the other 
hand clade ljs-bys-srs only has a posterior probability of 0.74 
(Fig. 3). An important feature of these trees was that the 
components of each clade showed a strong geographical associ- 
ation. All haplotypes of clade hby-lfy were from western localities, 
all haplotypes of clade nsc were from locality nsc, all haplotypes of 
clade ljs-bys-srs were from the eastern localities (Fig. 1). 

Population structure and genetic differentiation 

AMOVA analysis suggested that majority of the variation 
occurred among populations (82.15%) while differentiation within 
populations only contributed 17.85% to the total (Table 3). The 
fixation index, a measure of population differentiation due to 
genetic structure, was highly indicating a significant genetic 
differentiation among populations (/j-value =0.0±0.0) (Table 3). 



Table 4. F ST values between populations. 





Population 


bys 


hby 


Ify 




nsc srs 


bys 


0.0 










hby 


0.82921 (p = 0.0±0.0) 


0.0 








Ify 


0.89042 (p = 0.0±0.0) 


0.06828 (p = 0.0685 1 ± 0.002) 


0.0 






Ijs 


0.38150 (p = 0.0±0.0) 


0.84979 (p = 0.0±0.0) 


0.90412 (p = 0.0±0.0) 


0.0 




nsc 


0.86254 (p = 0.0±0.0) 


0.77573 (p = 0.0±0.0) 


0.83417 (p = 0.0±0.0) 


0.87731 (p = 0.0±0.0) 


0.0 


srs 


0.31488 (p = 0.0±0.0) 


0.79976 (p = 0.0±0.0) 


0.86822 (p = 0.0±0.0) 


0.33801 (p = 0.0±0.0) 


0.83596 (p = 0.0±0.0) 0.0 



doi:1 0.1 371 /joumal.pone.01 00729.t004 



PLOS ONE | www.plosone.org 



June 2014 | Volume 9 | Issue 6 | e100729 



Molecular Ecology of Scutiger ningshanensis 



Table 5. Geographical distances among populations. 





Population 


bys 


hby 


Ify 




nsc 


srs 


bys 


0.0 












hby 


400.14 


0.0 










Ify 


368.7 


38.76 


0.0 








Ijs 


20.05 


381.29 


350.31 


0.0 






nsc 


307.33 


103.54 


65.98 


289.6 


0.0 




srs 


40.36 


438.94 


408.07 


57.77 


347.22 


0.0 



doi:1 0.1 371 /journal.pone.01 00729.t005 



Population comparisons showed a significant genetic differen- 
tiation (F ST ) between most local populations except the population 
comparison between hby and Ify (Table 4). 

Mantel test results showed a significant correlation between the 
pairwise calculated genetic distance and pairwise calculated 
geographical distance (Table 5) of the populations (regression 
coefficient =0.001286, correlation coefficient =0.77051, p 
( r rand— r obs) = 0.0 1 85<0.05), indicating the existence of isolation- 
by-distance pattern of genetic divergence for cyt b+COI sequence, 
which suggests that the distribution of genetic variation is due to 
geographical separation rather than natural selection. The Mantel 
test results provided the evidence for large-scale geographical 
population structure in this species. 

Results from SAMOVA analysis indicated that the highest Fct 
equals to 0.84285 (p (rand.value > obs.value) = 0.01662±0.0) 
when K— 3, showing that the most likely number of populations is 
three. 

Historical biogeography 

The observed value of s is 55 which doesn't fall within the model 
of single-refugium indicating that the fragmentation model of 



population divergence was rejected (Fig. 4), that is, the current 
three lineages were derived from multiple refugia. 

Demographic history 

The results of mismatch distribution showed that the /(-values 
for SSDs and Rs of the total population and all clades were larger 
than 0.05, indicating a stable population size in the past; 
moreover, mismatch distribution in the total population and the 
clade hby-lfy showed frequency distribution of pairwise difference 
with multiple peaks (Fig. 5), on the other hand, the /(-values for 
Tajima's D of total population, clade hby-lfy and clade nsc were 
larger than 0.05 also indicating a stable population size in the past. 
However, the /(-values for Fu's Fgs of all clades were smaller than 
0.01, and the /(-values for Tajima's D of clade ljs-bys-srs was 
smaller than 0.05, indicating all clades experienced a sudden 
population expansion (Table 6); on the other hand, mismatch 
distribution in clades nsc and ljs-bys-srs showed frequency 
distribution of pairwise difference with single peak (Fig. 5), 
indicating an sudden population expansion in the past. The two 
statistics, Tajima's D value and the Fu's F s , are sensitive to 
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bottleneck effects or population expansion which leads to the more 
negative values of Tajim's D and Fu's F$ [71—74]. 

The allelic frequency spectrum for the entire population and all 
clades revealed an excess of singleton mutations and doesn't fit the 
neutral model (Fig. 6). The excess of singleton mutations could be 
caused by the expansion [58]. 

Bayesian skyline plot (BSP) further estimated the demographic 
history. The effective sample size (ESS) for all parameters of the 
BSP was greater than 200, showing that the 20 million generations 
were sufficient to determine the demographic history for each 
examined lineage. All lineages and the total population experi- 
enced population expansion. The hby-lfy lineage, nsc lineage and 
the total population experienced quick population growth after 
about 40,000 years ago, while the ljs-bys-srs lineage experienced a 
slow population growth after about 15,000 years ago. Noticeably, 
the hby-lfy lineage and the total population experienced a recent 
sharp decline in population, the ljs-bys-srs lineage showed a more 
recent population decline, while the nsc lineage maintained 
basically constant population size after the population expansion 
(Kg. 7). 

Genetic distances between lingeages 

The lineage hby-lfy and the lineage nsc were highly divergent 
from each other with an uncorrected p-distances of 4.2%, the 
lineage hby-lfy and ljs-bys-srs were also highly divergent with an 
uncorrected p-distances of 4.37%. Similarly, there was an high 
divergence between the lineage nsc and ljs-bys-srs, with an 
uncorrected p-distances of 2.37%. These values were similar to 
those p-distances (3%) reported between two different frog species 
[75]. Therefore S. ningshanensis probably contains three cryptic 
species or at least three subspecies. 

Discussion 

Genetic diversity 

Our results showed high levels of genetic diversity in S. 
ningshanensis. High-level genetic diversity in narrowly endemic 
species could be associated with the factors such as the long species 
history, the high-level gene flow, having experienced a recent 
contraction in population, multiple founder events, the mainte- 
nance of genetic diversity in refugia, the hybridization, and the 
ability to survive in a range of different habitats. [76-78]. All these 
factors except hybridization might contribute to the high level of 
genetic diversity in S. ningshanensis. This species has probably 
undergone a long history through which many mutations were 
accumulated since their origin. High-level gene flow is beneficial to 
fix the mutations in population. Multiple founder events were 
another alternative explanation for the high-level genetic diversity 
in this species. There is a possibility that multiple founder events, 
which accumulated more and more mutations, occurred in the 
history of this species. The two congeners of S. ningshanensis are S. 
boulengeri and S. liupanensis. The nearest localities of these two 
species are 167 km or 361 km respectively away from the locality 
hby of S. ningshanensis, which makes it impossible for the occurring 
of hybridization among these three species [79,80]. Therefore, 
hybridization may not contribute to the high-level genetic diversity 
in S. ningshanensis. The distributing region of S. ningshanensis include 
a range of different habitats such as high mountains with different 
elevation, streams, and different vegetations. High-level diversity 
in habitats also contributed to the high-level genetic diversity. 

Population structure 

Significant population structure occured based on the statistics 
pairwise differentiation. Most pairwise F^t values are high and 
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Figure 6. Allele frequency spectrum indicated an excess of singleton mutations in the combined mtDNA cyt b and COI sequences. 

Numbers above the line represent the number of sites with singleton mutations. 
doi:10.1371/journal.pone.0100729.g006 



statistically significant. Thus, the Ningshan alpine toad appears to 
exhibit substantial population differentiation across the Tsinling 
Mountains. 

Mantel test results showed a significant correlation between the 
genetic distance and geographical distance of the populations, 
indicating the presence of IBD pattern of genetic divergence for 
cyt b+COI sequences, suggesting that the distribution of genetic 
variation is due to geographical separation rather than natural 
selection. The Mantel test results provided the clear evidence for 
large-scale geographical population structure in S. ningshanensis. It 
is not possible that a significant Mantel test provided the evidence 
for discontinuity in the distribution of genetic variation. It rather 
showed a continuous distribution of the variation due to 
individuals mating preferentially with individuals from nearby 
populations [81]. 

AMOVA results indicated that 82.15% genetic variation 
occurred among populations, while differentiation within the 
populations only made 17.85% contributions. The high genetic 
variation among populations affirmed the presence of phylogeo- 
graphic structure in S. ningshanensis. SAMOVA results and the 
phylogenetic relationship analysis further affirmed the existence of 
phylogeographic structure in this species. 

Amphibians have poor dispersal capabilities and are sensitive to 
fine-scale landscape structure, topographic and altitudinal varia- 
tion and climatic changes [82-88]. Many amphibians are 
philopatric to breeding sites [89] which restricts gene flow and 
leads to significant genetic differentiation among populations and 
lineages. 



Topography and Pleistocene climate changes drive population 
genetic differentiation forming genetic structure pattern [90-92]. 
East Asia including China has undergone a series of cooler-drier 
climate cycles in the last 15 million years [93]. Dramatic climatic 
changes have caused the extinction of many organisms [94] . S. 
ningshanensis that distributed in the areas with low elevation 
disappeared during the interglacial in the Quaternary since it is an 
alpine species. In addition, S. ningshanensis retreated to a few refugia 
during glacial period. At the end of the Dali glaciation, S. 
ningshanensis experienced a rapid population expansion which 
occurred about 40,000 years before present, and it is experiencing 
a population contraction now probably due to anthropogenic 
disturbance. Topographic features, climatic fluctuation and 
anthropogenic activity together led to the current patchy 
geographic pattern for S. ningshanensis. This geographical distribu- 
tion pattern also seen in other montane organisms [95,96]. 

Finally, there is the possibility that unsampled populations 
between nsc and ljs may be genetically intermediate among the 
three groups. 

Historical biogeography 

The uplifts of the Tsinling Mountains promoted genetic 
differentiation among lineages of S. ningshanensis. The Tibetan 
Plateau experienced three phases of rapid uplifts [97,98]. The 
uplifts of the Tibetan Plateau impacted the environment of the 
surrounding areas including the Tsinling Mountains [97]. The 
Tsinling Mountains experienced similar uplift process, and three 
phases of uplifts occurred in Pliocene, Early Pleistocene and 
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Figure 7. Demographic patterns of each clade and the total population as determined from the Bayesian skyline plot (BSP). The X- 

axis is in units of million years in the past and the Y-axis is A/ e *u (effective population size x mutation rate per site per generation). The median 
estimates for the Iog10 of the population size are shown as thick solid lines, and the 95% highest posterior density (HPD) limits are shown by the 
shaded areas. 
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Holocene respectively [99-102]. Coalescent simulations indicated 
multiple refugia in S. ning.shanen.sis though we could not determine 
the number and location of the refugia based on the statistical 
phylogeography analyses. Most species of the genus Scutiger inhabit 
in the eastern Tibetan Plateau [81], we guess that S. ningshanensis 
evolved when the Scutiger species dispersed eastward along the 
Tsinling Mountains. The hby-lfy lineage split apart at first, split 
between the lineages nsc and ljs-bys-srs occurred when this species 
continued dispersing from west to east on the Tsinling Mountains, 
while no further split occurred after the second split. 

Demographic history 

Most demographic analyses revealed a sharp population 
expansion in all lineages of S. ningshanensis, and the expansion 
began simultaneously 40,000 years before present, which corre- 
sponds to the end of the Dali glaciation. Retreat of glacier led to 
population expansion in species on the Tsinling Mountains. A 
noteworthy phenomena is that S. ningshanensis is experiencing a 
distinct recent population contraction (though the lineage nsc is 
exceptional) as shown by the BSP analyses (Fig. 7), there is a 
possibility that anthropogenic disturbance resulted in the contrac- 
tion in population size of this species. Multiple uplifts of the 
Tsinling Mountains and fluctuations in population size of S. 
ningshanensis associated with glacial-interglacial cycles led to 



increases or decreases in the levels of genetic variation and 
coalescence times [23,25,28]. 

Putative cryptic species 

S. ningshanensis probably contains at least two cryptic species or 
subspecies based on the p-distance analysis though there is not 
difference in morphology among these cryptic species or 
subspecies. The geographical distances between the local popula- 
tions are long enough for occurrence of gene flow break in S. 
ningshanensis since it is philopatric to breeding sites. In addition, 
high peaks and deep valleys also contributed to the break of gene 
flow between populations. Subsequently, poor level of gene flow 
led to speciation events. However, there are limitations in 
taxonomic consequences based on only one taxonomic discipline. 
We will further confirm the cryptic speciation using data from 
nuclear DNA and ecological niche modeling (ENM) in the future. 

Conservation and management implications 

All potential cryptic species should be considered for conserva- 
tion. Different conservation strategies should be accepted for 
different species, therefore, it is inappropriate to protect a cryptic 
species complex using a single conservation strategy [103]. It is 
indispensable to understand and quantify biodiversity so that we 
can better explain and at last carry out conservation [103]. The 
distribution area of this species is limited to a narrow zone along 
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the Tsinling Mountains, if this species is a cryptic species complex 
itself, the distribution area of each cryptic species is much smaller 
which increases the risk for extinction. Effective protection 
measures are to be provided and carried out. Conservation should 
be considered for every cryptic species. 

Management units (MUs) are used to make conservation 
strategy for threatened species [104,105]. The three lineages 
(hby-lfy, nsc, ljs-bys-srs) derived from our study could be used as 
three different MUs, each MU has its own unique genetic 
structure. Every MU should be considered when conservation 
policy is made. 

Conclusions 

Three lineages were detected. The phylogenetic relationship 
analyses showed significant phylogeographic structure. The 
population expansion or contraction and genetic differentiation 
between populations or lineages could be explained by topography 
and the repetitive uplifts of the Tsinling Mountains and the 
climatic cycles during the late Pliocene and Pleistocene. S. 
ningshanensis experienced a rapid population expansion about 
40,000 years before present. The current decline in population size 
was probably caused by anthropogenic disturbance. Current 
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populations of S. ningshanensis are from different refugia though the 
location of these refugia could not be determined in our study. 
Topography is very important in shaping genetic connectivity. 
Topography, climatic changes and repetitive population expan- 
sion/contraction together led to the high level of genetic variation 
in S. ningshanensis. At least two putative cryptic species were 
detected. A total of three MUs were determined, which must be 
considered when conservation policy is made in future. 
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